'''
Created on Jun 19, 2011

@author: oabalbin
'''

from collections import defaultdict

def read_regions_file(rfile,ofile,gene2report):
    '''
    '''
    ifile = open(rfile)
    ofile = open(ofile,'w')
    ifile.next()
    
    genelist=defaultdict(list)
    
    for l in ifile:
        f = l.strip('\n').split('\t')
        chr, estart, eend, name, cov = f[0], f[1], f[2], f[3], float(f[6])
        name=name.split('|')[1] #+'|'+name.split('|')[0]
        genelist[name].append(cov)
        
    genes = sorted(genelist.keys())
    
    header=["geneName","avg_cov","min_cov","max_cov"]
    ofile.write(",".join(header).replace(',','\t')+'\n')
    for g in genes:
        covs = genelist[g]
    #for g,covs in genelist.iteritems():
        average = float(sum(covs)) / len(covs)
        l =map(str,[g,average, int(min(covs)), int(max(covs))])
        if g in gene2report:
            ofile.write(",".join(l).replace(',','\t')+'\n')
        else:
            print l
           
def informative_genes_report(genes):
    '''
    chr2    30142808        30144482        NM_004304|ALK
    '''
    thgens=[]
    gfile=open(genes)
    
    for g in gfile:
        f = g.strip('\n').split('\t')
        name=f[3].split('|')[1]
        thgens.append(name)
        
    return thgens

genes = informative_genes_report('/exds/users/oabalbin/projects/exomes/informative_genes_annotation.bed')
read_regions_file('/exds/users/oabalbin/projects/exomes/tmp/MO1001_n_informative_genes_coverage.bed', 
                  '/exds/users/oabalbin/projects/exomes/tmp/MO1001_n_informative_genes_coverage.report',
                  genes)
